{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Refinery\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "In this example, we’ll demonstrate how you can use mathematical optimization to optimize the output of a refinery. You’ll learn how to generate an optimal production plan that maximizes total profit, while taking into account production capacity and other restrictions.\n",
    "\n",
    "More information on this type of model can be found in example # 6 of the fifth edition of Modeling Building in Mathematical Programming by H. P. Williams on pages 258 and 306 – 310.\n",
    "\n",
    "This modeling example is at the intermediate level, where we assume that you know Python and are familiar with the Gurobi Python API. In addition, you should have some knowledge about building mathematical optimization models.\n",
    "\n",
    "**Download the Repository** <br />\n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). \n",
    "\n",
    "---\n",
    "## Problem Description\n",
    "\n",
    "An oil refinery purchases two types of crude oil (Crude 1 and Crude 2) and refines them through a four-step process (distillation, reforming, cracking, and blending) to produce gas and fuels that are sold.\n",
    "\n",
    "### Step One: Distillation\n",
    "\n",
    "The distillation process separates each crude oil into six fractions according to their boiling points. These fractions are: light naphtha, medium naphtha, heavy naphtha, light oil, heavy oil, and residuum. The octane numbers for the light, medium and heavy naphthas are, respectively, 90, 80, and 70.\n",
    "\n",
    "One barrel of crude separates into following fractions:\n",
    "\n",
    "| <i></i> | Light Naphta | Medium Naphta | Heavy Naphta | Light Oil | Heavy Oil | Residuum |\n",
    "| --- | --- | --- | --- | --- | --- | --- |\n",
    "| Crude 1 | 0.1 | 0.2 | 0.2 | 0.12 | 0.2 | 0.13 |\n",
    "| Crude 2 | 0.15 | 0.25 | 0.18 | 0.08 | 0.19 | 0.12 |\n",
    "\n",
    "There is a small amount of waste in distillation (5% and 3% for Crude 1 and 2, respectively).\n",
    "\n",
    "### Step Two: Reforming\n",
    "\n",
    "After distillation, the resulting naphthas can be blended into different grades of gas, or they can go through a process called reforming. The output of the reforming process is a product known as reformed gasoline with an octane number of 115.\n",
    "\n",
    "The yields of reformed gasoline from each barrel of the different naphthas are given as follows:\n",
    "\n",
    "- One barrel of light naphtha yields 0.6 barrels of reformed gasoline.\n",
    "- One barrel of medium naphtha yields 0.52 barrels of reformed gasoline.\n",
    "- One barrel of heavy naphtha yields 0.45 barrels of reformed gasoline.\n",
    "\n",
    "### Step Three: Cracking\n",
    "\n",
    "Light and heavy oils can be blended into jet fuel or be put through a process known as catalytic cracking. The catalytic cracker produces cracked oil and cracked gasoline.\n",
    "\n",
    "Cracked gasoline has an octane number of 105 with the following yields:\n",
    "\n",
    "- One barrel of light oil yields 0.68 barrels of cracked oil and 0.28 barrels of cracked gasoline.\n",
    "- One barrel of heavy oil yields 0.75 barrels of cracked oil and 0.2 barrels of cracked gasoline.\n",
    "\n",
    "Cracked oil is used for blending fuel oil and jet fuel; cracked gasoline is used for blending gasoline. Residuum can be used for either producing lube-oil or blending into jet fuel and fuel oil:\n",
    "\n",
    "- One barrel of residuum yields 0.5 barrels of lube-oil.\n",
    "\n",
    "### Step Four: Blending\n",
    "\n",
    "#### Gasoline\n",
    "\n",
    "There are two kinds of gasoline — regular and premium — made by blending naphtha, reformed gasoline, and cracked gasoline. The only requirement is that regular gasoline must have an octane of at least 84 and premium gasoline must have an octane number of at least 94. It is assumed that octane numbers blend linearly by volume.\n",
    "\n",
    "#### Jet Fuel\n",
    "\n",
    "Jet fuel must have a vapor pressure that does not exceed $1.0 \\frac{\\text{kg}}{\\text{cm}^2}$. The vapor pressures for light, heavy, cracked oils and residuum are, respectively, 1.0, 0.6, 1.5 and 0.05 $\\frac{\\text{kg}}{\\text{cm}^2}$. It is assumed that vapor pressures blend linearly by volume.\n",
    "\n",
    "#### Fuel Oil\n",
    "\n",
    "To produce fuel oil, you must blend light oil, cracked oil, heavy oil and residuum in the ratio of 10:4:3:1. \n",
    "\n",
    "The availability and capacity limitations are as follows:\n",
    "\n",
    "- The daily availability of crude 1 is 20 000 barrels.\n",
    "- The daily availability of crude 2 is 30 000 barrels.\n",
    "- At most 45 000 barrels of crude can be distilled per day.\n",
    "- At most 10 000 barrels of naphtha can be reformed per day.\n",
    "- At most 8000 barrels of oil can be cracked per day.\n",
    "- The daily production of lube oil must be between 500 and 1000 barrels.\n",
    "- Premium gasoline production must be at least 40% of regular gasoline production.\n",
    "\n",
    "Each final product has the following profit contribution per barrel (pennies per barrel):\n",
    "\n",
    "| <i></i> | Profit Contribution |\n",
    "| --- | --- |\n",
    "| Premium Gasoline | 700 |\n",
    "| Regular Gasoline | 600 |\n",
    "| Jet Fuel | 400 |\n",
    "| Fuel Oil | 350 |\n",
    "| Lube Oil | 150 |\n",
    "\n",
    "The key question is: How should the operations of the refinery be planned in order to maximize total profit?\n",
    "\n",
    "---\n",
    "## Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "\n",
    "$i \\in \\text{Crudes}=\\{1,2\\}$: Set of crude oils.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$\\text{buy_limit}_i \\in \\mathbb{N}$: Maximum number of barrels of crude $i$ to buy.\n",
    "\n",
    "$\\text{distill_cap} \\in \\mathbb{N}$: Maximum number of barrels of crude oil to distill. \n",
    "\n",
    "$\\text{reform_cap} \\in \\mathbb{N}$: Maximum number of barrels of naphtha to reform. \n",
    "\n",
    "$\\text{crack_cap} \\in \\mathbb{N}$ Maximum number of barrels of oil to crack.\n",
    "\n",
    "$\\text{LBO_min}, \\text{LBO_max} \\in \\mathbb{N}$ Minimum and Maximum number of barrels of lube oil to produce.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$\\text{CR}_i \\in [0,\\text{buy_limit}_i] \\subset \\mathbb{R}^+$: Number of barrels of crude $i$ to buy.\n",
    "\n",
    "$\\text{LN} \\in \\mathbb{R}^+$: Number of barrels of light naphtha to distill.\n",
    "\n",
    "$\\text{MN} \\in \\mathbb{R}^+$: Number of barrels of medium naphtha to distill.\n",
    "\n",
    "$\\text{HN} \\in \\mathbb{R}^+$: Number of barrels of heavy naphtha to distill.\n",
    "\n",
    "$\\text{LO} \\in \\mathbb{R}^+$: Number of barrels of light oil to distill.\n",
    "\n",
    "$\\text{HO} \\in \\mathbb{R}^+$: Number of barrels of heavy oil to distill.\n",
    "\n",
    "$\\text{R} \\in \\mathbb{R}^+$: Number of barrels of residuum to distill.\n",
    "\n",
    "$\\text{LNRG} \\in \\mathbb{R}^+$: Number of barrels of light naphtha used to produce reformed gasoline.\n",
    "\n",
    "$\\text{MNRG} \\in \\mathbb{R}^+$: Number of barrels of medium naphtha used to produce reformed gasoline.\n",
    "\n",
    "$\\text{HNRG} \\in \\mathbb{R}^+$: Number of barrels of heavy naphtha used to produce reformed gasoline.\n",
    "\n",
    "$\\text{RG} \\in \\mathbb{R}^+$: Number of barrels of reformed gasoline to produce.\n",
    "\n",
    "$\\text{LOCGO} \\in \\mathbb{R}^+$: Number of barrels of light oil used to produce cracked gasoline and cracked oil.\n",
    "\n",
    "$\\text{HOCGO} \\in \\mathbb{R}^+$: Number of barrels of heavy oil used to produce cracked gasoline and cracked oil.\n",
    "\n",
    "$\\text{CG} \\in \\mathbb{R}^+$: Number of barrels of cracked gasoline to produce.\n",
    "\n",
    "$\\text{CO} \\in \\mathbb{R}^+$: Number of barrels of cracked oil to produce.\n",
    "\n",
    "$\\text{LNPMF} \\in \\mathbb{R}^+$: Number of barrels of light naphtha used to produce premium motor fuel.\n",
    "\n",
    "$\\text{LNRMF} \\in \\mathbb{R}^+$: Number of barrels of light naphtha used to produce regular motor fuel.\n",
    "\n",
    "$\\text{MNPMF} \\in \\mathbb{R}^+$: Number of barrels of medium naphtha used to produce premium motor fuel.\n",
    "\n",
    "$\\text{MNRMF} \\in \\mathbb{R}^+$: Number of barrels of medium naphtha used to produce regular motor fuel.\n",
    "\n",
    "$\\text{HNPMF} \\in \\mathbb{R}^+$: Number of barrels of heavy naphtha used to produce premium motor fuel.\n",
    "\n",
    "$\\text{HNRMF} \\in \\mathbb{R}^+$: Number of barrels of heavy naphtha used to produce regular motor fuel.\n",
    "\n",
    "$\\text{RGPMF} \\in \\mathbb{R}^+$: Number of barrels of reformed gasoline used to produce premium motor fuel.\n",
    "\n",
    "$\\text{RGRMF} \\in \\mathbb{R}^+$: Number of barrels of reformed gasoline used to produce regular motor fuel.\n",
    "\n",
    "$\\text{CGPMF} \\in \\mathbb{R}^+$: Number of barrels of cracked gasoline used to produce premium motor fuel.\n",
    "\n",
    "$\\text{CGRMF} \\in \\mathbb{R}^+$: Number of barrels of cracked gasoline used to produce regular motor fuel.\n",
    "\n",
    "$\\text{LOJF} \\in \\mathbb{R}^+$: Number of barrels of light oil used to produce jet fuel.\n",
    "\n",
    "$\\text{HOJF} \\in \\mathbb{R}^+$: Number of barrels of heavy oil used to produce jet fuel.\n",
    "\n",
    "$\\text{RJF} \\in \\mathbb{R}^+$: Number of barrels of residuum used to produce jet fuel.\n",
    "\n",
    "$\\text{COJF} \\in \\mathbb{R}^+$: Number of barrels of cracked oil used to produce jet fuel.\n",
    "\n",
    "$\\text{RLBO} \\in \\mathbb{R}^+$: Number of barrels of residuum used to produce lube oil.\n",
    "\n",
    "$\\text{PMF} \\in \\mathbb{R}^+$: Number of barrels of premium motor fuel to produce.\n",
    "\n",
    "$\\text{RMF} \\in \\mathbb{R}^+$: Number of barrels of regular motor fuel to produce.\n",
    "\n",
    "$\\text{JF} \\in \\mathbb{R}^+$: Number of barrels of jet fuel to produce.\n",
    "\n",
    "$\\text{FO} \\in \\mathbb{R}^+$: Number of barrels of fuel oil to produce.\n",
    "\n",
    "$\\text{LBO} \\in [\\text{LBO_min}, \\text{LBO_max}] \\subset \\mathbb{R}^+$: Number of barrels of lube oil to produce.\n",
    "\n",
    "### Objective Function\n",
    "\n",
    "- **Profit:** Maximize the total profit (in hundreds of USD).\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Max} \\quad Z = 7*\\text{PMF} + 6*\\text{RMF} + 4*\\text{JF} + 3.5*\\text{FO} + 1.5*\\text{LBO}\n",
    "\\tag{0}\n",
    "\\end{equation}\n",
    "\n",
    "### Constraints\n",
    "\n",
    "- **Distillation Capacity**: The number of barrels of crude oil to distill cannot exceed the capacity.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in \\text{Crudes}}{\\text{CR}_i} \\leq \\text{distill_cap}\n",
    "\\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- **Reforming Capacity**: The number of barrels of naphtha to reform cannot exceed the capacity.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LNRG} + \\text{MNRG} + \\text{HNRG} \\leq \\text{reform_cap}\n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "- **Cracking Capacity**: The number of barrels of oil to crack cannot exceed the capacity.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LOCGO} + \\text{HOCGO} \\leq \\text{crack_cap}\n",
    "\\tag{3}\n",
    "\\end{equation}\n",
    "\n",
    "- **Yield**: The number of barrels produced depends on the quantities of inputs used, as well as their corresponding yields.\n",
    "\n",
    "\\begin{equation}\n",
    "0.10*\\text{CR}_1 + 0.15*\\text{CR}_2 = \\text{LN}\n",
    "\\tag{4.1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.20*\\text{CR}_1 + 0.25*\\text{CR}_2 = \\text{MN}\n",
    "\\tag{4.2}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.20*\\text{CR}_1 + 0.18*\\text{CR}_2 = \\text{HN}\n",
    "\\tag{4.3}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.12*\\text{CR}_1 + 0.08*\\text{CR}_2 = \\text{LO}\n",
    "\\tag{4.4}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.20*\\text{CR}_1 + 0.19*\\text{CR}_2 = \\text{HO}\n",
    "\\tag{4.5}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.13*\\text{CR}_1 + 0.12*\\text{CR}_2 = \\text{R}\n",
    "\\tag{4.6}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.60*\\text{LNRG} + 0.52*\\text{MNRG} + 0.45*\\text{HNRG} = \\text{RG}\n",
    "\\tag{4.7}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.68*\\text{LOCGO} + 0.75*\\text{HOCGO} = \\text{CO}\n",
    "\\tag{4.8}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.28*\\text{LOCGO} + 0.20*\\text{HOCGO} = \\text{CG}\n",
    "\\tag{4.9}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "0.50*\\text{RLBO} = \\text{LBO}\n",
    "\\tag{4.10}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LNPMF} + \\text{MNPMF} + \\text{HNPMF} + \\text{RGPMF} + \\text{CGPMF} = \\text{PMF}\n",
    "\\tag{4.11}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LNRMF} + \\text{MNRMF} + \\text{HNRMF} + \\text{RGRMF} + \\text{CGRMF} = \\text{RMF}\n",
    "\\tag{4.12}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LOJF} + \\text{HOJF} + \\text{COJF} + \\text{RJF} = \\text{JF}\n",
    "\\tag{4.13}\n",
    "\\end{equation}\n",
    "\n",
    "- **Mass Conservation:** The number of barrels used must be equal to the number of barrels available.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LNRG} + \\text{LNPMF} + \\text{LNRMF} = \\text{LN}\n",
    "\\tag{5.1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{MNRG} + \\text{MNPMF} + \\text{MNRMF} = \\text{MN}\n",
    "\\tag{5.2}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{HNRG} + \\text{HNPMF} + \\text{HNRMF} = \\text{HN}\n",
    "\\tag{5.3}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{LOCGO} + \\text{LOJF} + 0.55*\\text{FO} = \\text{LO}\n",
    "\\tag{5.4}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{HOCGO} + \\text{HOJF} + 0.17*\\text{FO} = \\text{HO}\n",
    "\\tag{5.5}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{COJF} + 0.22*\\text{FO} = \\text{CO}\n",
    "\\tag{5.6}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{RLBO} + \\text{RJF} + 0.0555*\\text{FO} = \\text{R}\n",
    "\\tag{5.7}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{CGPMF} + \\text{CGRMF} = \\text{CG}\n",
    "\\tag{5.8}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{RGPMF} + \\text{RGRMF} = \\text{RG}\n",
    "\\tag{5.9}\n",
    "\\end{equation}\n",
    "\n",
    "- **Premium-to-Regular Proportion:** The production ratio between premium and regular gasoline must satisfy the minimum requirement.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{PMF} \\geq 0.40*\\text{RMF}\n",
    "\\tag{6}\n",
    "\\end{equation}\n",
    "\n",
    "- **Octane Tolerance:** The octane rating of each type of gasoline cannot drop below the lower tolerance.\n",
    "\n",
    "\\begin{equation}\n",
    "90*\\text{LNPMF} + 80*\\text{MNPMF} + 70*\\text{HNPMF} + 115*\\text{RGPMF} + 105*\\text{CGPMF} \\geq 94*\\text{PMF}\n",
    "\\tag{7.1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "90*\\text{LNRMF} + 80*\\text{MNRMF} + 70*\\text{HNRMF} + 115*\\text{RGRMF} + 105*\\text{CGRMF} \\geq 84*\\text{PMF}\n",
    "\\tag{7.2}\n",
    "\\end{equation}\n",
    "\n",
    "- **Vapor-Pressure Tolerance:** The vapor pressure of jet fuel cannot drop below the lower tolerance.\n",
    "\n",
    "\\begin{equation}\n",
    "1.0*\\text{LOJF} + 0.6*\\text{HOJF} + 1.5*\\text{COJF} + 0.05*\\text{RJF} \\leq 1.0*\\text{JF}\n",
    "\\tag{8}\n",
    "\\end{equation}\n",
    "\n",
    "---\n",
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.11 & Gurobi 11.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input Data\n",
    "We define all the input data of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "\n",
    "crude_numbers = range(1,2+1)\n",
    "petrols = [\"Premium_fuel\", \"Regular_fuel\"]\n",
    "end_product_names = [\"Premium_fuel\", \"Regular_fuel\", \"Jet_fuel\", \"Fuel_oil\", \"Lube_oil\"]\n",
    "distillation_products_names = [\"Light_naphtha\", \"Medium_naphtha\", \"Heavy_naphtha\",\n",
    "                               \"Light_oil\", \"Heavy_oil\", \"Residuum\"]\n",
    "naphthas = [\"Light_naphtha\", \"Medium_naphtha\", \"Heavy_naphtha\"]\n",
    "intermediate_oils = [\"Light_oil\", \"Heavy_oil\"]\n",
    "cracking_products_names = [\"Cracked_gasoline\", \"Cracked_oil\"]\n",
    "used_for_motor_fuel_names = [\"Light_naphtha\", \"Medium_naphtha\", \"Heavy_naphtha\",\n",
    "                             \"Reformed_gasoline\", \"Cracked_gasoline\"]\n",
    "used_for_jet_fuel_names = [\"Light_oil\", \"Heavy_oil\", \"Residuum\", \"Cracked_oil\"]\n",
    "\n",
    "buy_limit = {1:20000, 2:30000}\n",
    "lbo_min = 500\n",
    "lbo_max = 1000\n",
    "\n",
    "distill_cap = 45000\n",
    "reform_cap = 10000\n",
    "crack_cap = 8000\n",
    "\n",
    "distillation_splitting_coefficients = {\"Light_naphtha\": (0.1, 0.15),\n",
    "                          \"Medium_naphtha\": (0.2, 0.25),\n",
    "                         \"Heavy_naphtha\": (0.2, 0.18),\n",
    "                         \"Light_oil\": (0.12, 0.08),\n",
    "                         \"Heavy_oil\": (0.2, 0.19),\n",
    "                         \"Residuum\": (0.13, 0.12)}\n",
    "\n",
    "cracking_splitting_coefficients = {(\"Light_oil\",\"Cracked_oil\"): 0.68,\n",
    "                                   (\"Heavy_oil\",\"Cracked_oil\"): 0.75,\n",
    "                                   (\"Light_oil\",\"Cracked_gasoline\"): 0.28,\n",
    "                                   (\"Heavy_oil\",\"Cracked_gasoline\"): 0.2}\n",
    "\n",
    "reforming_splitting_coefficients = {\"Light_naphtha\": 0.6, \"Medium_naphtha\":0.52, \"Heavy_naphtha\":0.45}\n",
    "end_product_profit = {\"Premium_fuel\":7, \"Regular_fuel\":6, \"Jet_fuel\":4, \"Fuel_oil\":3.5, \"Lube_oil\":1.5}\n",
    "blending_coefficients = {\"Light_oil\": 0.55, \"Heavy_oil\": 0.17, \"Cracked_oil\": 0.22, \"Residuum\": 0.055}\n",
    "\n",
    "lube_oil_factor = 0.5\n",
    "pmf_rmf_ratio = 0.4\n",
    "\n",
    "octance_number_coefficients = {\n",
    "    \"Light_naphtha\":90,\n",
    "    \"Medium_naphtha\":80,\n",
    "    \"Heavy_naphtha\":70,\n",
    "    \"Reformed_gasoline\":115,\n",
    "    \"Cracked_gasoline\":105,\n",
    "}\n",
    "octance_number_fuel = {\"Premium_fuel\": 94,\"Regular_fuel\": 84}\n",
    "\n",
    "vapor_pressure_constants = [0.6, 1.5, 0.05]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "We create a model and the variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n"
     ]
    }
   ],
   "source": [
    "refinery = gp.Model('Refinery_Optimization')\n",
    "\n",
    "# Variables\n",
    "crudes = refinery.addVars(crude_numbers, ub=buy_limit, name=\"cr\")    \n",
    "end_products = refinery.addVars(end_product_names, name=\"end_prod\")\n",
    "end_products[\"Lube_oil\"].lb= lbo_min\n",
    "end_products[\"Lube_oil\"].ub= lbo_max\n",
    "distillation_products = refinery.addVars(distillation_products_names, name=\"dist_prod\")\n",
    "reform_usage = refinery.addVars(naphthas, name=\"napthas_to_reformed_gasoline\")\n",
    "reformed_gasoline = refinery.addVar(name=\"reformed_gasoline\")\n",
    "cracking_usage = refinery.addVars(intermediate_oils,name=\"intermediate_oils_to_cracked_gasoline\")\n",
    "cracking_products = refinery.addVars(cracking_products_names,  name=\"cracking_prods\")\n",
    "used_for_regular_motor_fuel = refinery.addVars(used_for_motor_fuel_names, name=\"motor_fuel_to_regular_motor_fuel\")\n",
    "used_for_premium_motor_fuel = refinery.addVars(used_for_motor_fuel_names, name=\"motot_fuel_to_premium_motor_fuel\")\n",
    "used_for_jet_fuel = refinery.addVars(used_for_jet_fuel_names, name=\"jet_fuel\")\n",
    "used_for_lube_oil = refinery.addVar(vtype=GRB.CONTINUOUS,name=\"residuum_used_for_lube_oil\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we insert the constraints.\n",
    "\n",
    "The distillation capacity constraint is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Distillation capacity\n",
    "DistillationCap = refinery.addConstr(crudes.sum() <= distill_cap, \"Distill_cap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reforming capacity constraint is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Reforming capacity\n",
    "ReformingCap = refinery.addConstr(reform_usage.sum() <= reform_cap, \"Reform_cap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cracking capacity constraint is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3. Cracking capacity\n",
    "CrackingCap = refinery.addConstr(cracking_usage.sum() <= crack_cap, \"Crack_cap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantity of distillation products produced depends on the quantity of crude oil used, taking into account the way in which each crude splits under distillation. This gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.1-4.6 Yield (Crude oil products)\n",
    "YieldCrudeOil = refinery.addConstrs((gp.quicksum(distillation_splitting_coefficients[dpn][crude-1]*crudes[crude] for crude in crudes)\n",
    "                  == distillation_products[dpn] for dpn in distillation_products_names), \"Splitting_distillation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantity of reformed gasoline produced depends on the quantities of naphthas used in the reforming process. This gives the constraint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.7 Yield (Reforming of Naphthas)\n",
    "YieldNaphthas = refinery.addConstr(reform_usage.prod(reforming_splitting_coefficients) == reformed_gasoline, \"Splitting_reforming\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantities of cracked oil and cracked gasoline produced depend on the quantities of light and heavy oil used. This gives the constraints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.8-4.9 Yield (Cracking of oils)\n",
    "YieldCrackingOil = refinery.addConstrs((gp.quicksum(cracking_splitting_coefficients[oil, crack_prod]*cracking_usage[oil]\n",
    "                           for oil in intermediate_oils) == cracking_products[crack_prod]\n",
    "                  for crack_prod in cracking_products_names),\n",
    "                 name=\"Splitting_cracking\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantity of lube-oil produced (and sold) is 0.5 times the quantity of residuum used. This gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.10 Yield (Lube oil)\n",
    "YieldLubeOil = refinery.addConstr(lube_oil_factor*used_for_lube_oil == end_products[\"Lube_oil\"],\n",
    "                \"Splitting_lube_oil\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantity of motor fuels and jet fuel produced is equal to the total quantity of their ingredients. This gives the constraints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4.11 Yield (Premium gasoline)\n",
    "YieldPremium = refinery.addConstr(used_for_premium_motor_fuel.sum() == end_products[\"Premium_fuel\"], \"Blending_premium_fuel\")\n",
    "\n",
    "# 4.12 Yield (Regular gasoline)\n",
    "YieldRegular = refinery.addConstr(used_for_regular_motor_fuel.sum() == end_products[\"Regular_fuel\"], \"Blending_regular_fuel\")\n",
    "\n",
    "# 4.13 Yield (Jet fuel)\n",
    "YieldJetFuel = refinery.addConstr(used_for_jet_fuel.sum() == end_products[\"Jet_fuel\"], \"Continuity_jet_fuel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantities of naphthas used for reforming and blending are equal to the quantities available. This gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5.1-5.3 Mass conservation (Naphthas)\n",
    "MassBalNaphthas = refinery.addConstrs((reform_usage[naphtha] +\n",
    "                    used_for_regular_motor_fuel[naphtha] +\n",
    "                    used_for_premium_motor_fuel[naphtha] ==\n",
    "                    distillation_products[naphtha] for naphtha in naphthas), \"Continuity_napththa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the blending of fuel oil, the proportion of light oil/heavy oil/ cracked oil/ residuum is fixed. Therefore, separate variables have not been introduced for this proportion as it is determined by the variables. This gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5.4 Mass Conservation (Light oil)\n",
    "MassBalLightOil = refinery.addConstr(cracking_usage[\"Light_oil\"]+\n",
    "                used_for_jet_fuel[\"Light_oil\"]+\n",
    "                blending_coefficients[\"Light_oil\"]*end_products[\"Fuel_oil\"] ==\n",
    "                distillation_products[\"Light_oil\"], \"Fixed_proportion_light_oil_for_blending\")\n",
    "\n",
    "# 5.5 Mass Conservation (Heavy oil)\n",
    "MassBalHeavyOil = refinery.addConstr(cracking_usage[\"Heavy_oil\"]+\n",
    "                used_for_jet_fuel[\"Heavy_oil\"]+\n",
    "                blending_coefficients[\"Heavy_oil\"]*end_products[\"Fuel_oil\"] ==\n",
    "                distillation_products[\"Heavy_oil\"], \"Fixed_proportion_heavy_oil_for_blending\")\n",
    "\n",
    "# 5.6 Mass Conservation (Cracked oil)\n",
    "MassBalCrackedOil = refinery.addConstr(used_for_jet_fuel[\"Cracked_oil\"]+\n",
    "                blending_coefficients[\"Cracked_oil\"]*end_products[\"Fuel_oil\"] ==\n",
    "                cracking_products[\"Cracked_oil\"], \"Fixed_proportion_cracked_oil_for_blending\")\n",
    "\n",
    "# 5.7 Mass Conservation (Residuum)\n",
    "MassBalResiduum = refinery.addConstr(used_for_lube_oil +\n",
    "                used_for_jet_fuel[\"Residuum\"]+\n",
    "                blending_coefficients[\"Residuum\"]*end_products[\"Fuel_oil\"] ==\n",
    "                distillation_products[\"Residuum\"], \"Fixed_proportion_residuum_for_blending\")\n",
    "\n",
    "# 5.8 Mass conservation (Cracked gasoline)\n",
    "MassBalCrackedGas = refinery.addConstr(used_for_regular_motor_fuel[\"Cracked_gasoline\"] +\n",
    "                used_for_premium_motor_fuel[\"Cracked_gasoline\"] ==\n",
    "                cracking_products[\"Cracked_gasoline\"], \"Continuity_cracked_gasoline\")\n",
    "\n",
    "# 5.9 Mass conservation (Reformed gasoline)\n",
    "MassBalReformedGas = refinery.addConstr(used_for_regular_motor_fuel[\"Reformed_gasoline\"] +\n",
    "                used_for_premium_motor_fuel[\"Reformed_gasoline\"] ==\n",
    "                reformed_gasoline, \"Continuity_reformed_gasoline\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Premium motor fuel production must be at least 40% of regular motor fuel production, giving:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 7. Premium-to-regular proportion\n",
    "Premium2Regular = refinery.addConstr(end_products[\"Premium_fuel\"] >= pmf_rmf_ratio*end_products[\"Regular_fuel\"],\n",
    "                \"Prem2reg_prop\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is necessary to stipulate that the octane number of premium motor (regular motor) fuel does not drop below 94 (84). This is done by the constraint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 8.1-8.2 Octane tolerance\n",
    "OctaneRegular = refinery.addConstr(used_for_regular_motor_fuel.prod(octance_number_coefficients) >=\n",
    "                octance_number_fuel[\"Regular_fuel\"] * end_products[\"Regular_fuel\"],\n",
    "                \"Octane_tol_regular_fuel\")\n",
    "\n",
    "OctanePremium = refinery.addConstr(used_for_premium_motor_fuel.prod(octance_number_coefficients) >=\n",
    "                octance_number_fuel[\"Premium_fuel\"] * end_products[\"Premium_fuel\"],\n",
    "                \"Octane_tol_premium_fuel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For jet fuel, we have the constraint imposed by vapor pressure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9. Vapor-pressure tolerance\n",
    "VaporPressure = refinery.addConstr(used_for_jet_fuel[\"Light_oil\"] +\n",
    "                vapor_pressure_constants[0]*used_for_jet_fuel[\"Heavy_oil\"] +\n",
    "                vapor_pressure_constants[1]*used_for_jet_fuel[\"Cracked_oil\"] +\n",
    "                vapor_pressure_constants[2]*used_for_jet_fuel[\"Residuum\"] <= end_products[\"Jet_fuel\"],\n",
    "                \"Vapor_pressure_tol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This model has 29 constraints together with simple bounds on three variables.\n",
    "\n",
    "A note should be made concerning the blending of fuel oil where the ingredients (light, heavy, and cracked oil and residuum) are used in fixed proportions. It might be preferable to think of the production of fuel oil as an activity. It is common in the oil industry to think in terms of activities rather than quantities. In Section 3.4 of the H.P. Williams book, model formulations are discussed where activities represent the extreme modes of operation of a process. In this example, we have a special case of a process with one mode of operation. The level of this activity then fixes the ratios of the ingredients, in a case such as this, automatically.\n",
    "\n",
    "The only variables involving a profit (or cost) are the final products. This gives an objective (in dollars) to be maximized of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 0. Profit\n",
    "refinery.setObjective(end_products.prod(end_product_profit), GRB.MAXIMIZE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the optimization process starts and Gurobi finds the optimal solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 29 rows, 36 columns and 106 nonzeros\n",
      "Model fingerprint: 0xe30699e6\n",
      "Coefficient statistics:\n",
      "  Matrix range     [5e-02, 1e+02]\n",
      "  Objective range  [2e+00, 7e+00]\n",
      "  Bounds range     [5e+02, 3e+04]\n",
      "  RHS range        [8e+03, 5e+04]\n",
      "Presolve removed 13 rows and 14 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 16 rows, 22 columns, 72 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    1.1887574e+06   6.045565e+04   0.000000e+00      0s\n",
      "      14    2.1136513e+05   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 14 iterations and 0.01 seconds\n",
      "Optimal objective  2.113651348e+05\n"
     ]
    }
   ],
   "source": [
    "refinery.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Analysis\n",
    "\n",
    "The optimal solution results in a profit of $\\$211,365.13$. The optimal values of the variables are given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cr[1] = 15000.0\n",
      "cr[2] = 30000.0\n",
      "end_prod[Premium_fuel] = 6817.78\n",
      "end_prod[Regular_fuel] = 17044.45\n",
      "end_prod[Jet_fuel] = 15156.0\n",
      "end_prod[Lube_oil] = 500.0\n",
      "dist_prod[Light_naphtha] = 6000.0\n",
      "dist_prod[Medium_naphtha] = 10500.0\n",
      "dist_prod[Heavy_naphtha] = 8400.0\n",
      "dist_prod[Light_oil] = 4200.0\n",
      "dist_prod[Heavy_oil] = 8700.0\n",
      "dist_prod[Residuum] = 5550.0\n",
      "napthas_to_reformed_gasoline[Heavy_naphtha] = 5406.86\n",
      "reformed_gasoline = 2433.09\n",
      "intermediate_oils_to_cracked_gasoline[Light_oil] = 4200.0\n",
      "intermediate_oils_to_cracked_gasoline[Heavy_oil] = 3800.0\n",
      "cracking_prods[Cracked_gasoline] = 1936.0\n",
      "cracking_prods[Cracked_oil] = 5706.0\n",
      "motor_fuel_to_regular_motor_fuel[Light_naphtha] = 273.07\n",
      "motor_fuel_to_regular_motor_fuel[Medium_naphtha] = 10500.0\n",
      "motor_fuel_to_regular_motor_fuel[Heavy_naphtha] = 2993.14\n",
      "motor_fuel_to_regular_motor_fuel[Reformed_gasoline] = 1342.24\n",
      "motor_fuel_to_regular_motor_fuel[Cracked_gasoline] = 1936.0\n",
      "motot_fuel_to_premium_motor_fuel[Light_naphtha] = 5726.93\n",
      "motot_fuel_to_premium_motor_fuel[Reformed_gasoline] = 1090.84\n",
      "jet_fuel[Heavy_oil] = 4900.0\n",
      "jet_fuel[Residuum] = 4550.0\n",
      "jet_fuel[Cracked_oil] = 5706.0\n",
      "residuum_used_for_lube_oil = 1000.0\n"
     ]
    }
   ],
   "source": [
    "for var in refinery.getVars():\n",
    "    if abs(var.x) > 1e-6:\n",
    "        print(\"{0} = {1}\".format(var.varName, np.round(var.x, 2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
